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Abstract 

A sensitivity analysis of general stoichiometric networks is considered. The results are 
presented as a generalization of Metabolic Control Analysis, which has been concerned pri- 
marily with system sensitivities at steady state. An expression for time-varying sensitivity 
coefficients is given, and the Summation and Connectivity Theorems are generalized. The 
results are compared to previous treatments. The analysis is accompanied by a discussion 
of the computation of the sensitivity coefficients and an application to a model of photo- 
transduction. 



1 Introduction 

Sensitivity analysis is an important tool in the study of systems dependent on external param- 
eters. By their nature, stoichiometric networks allow an elegant description of the relationships 
between component and systemic sensitivities. These relationships, particularly those existing 
at equilibria of the system, have been addressed by the fields of Biochemical Systems Theory 
(BST) and Metabolic Control Analysis (MCA) g §. 



This paper describes the computation and interpretation of system sensitivities over arbi- 
trary trajectories. The results are presented in the framework provided by MCA, which is a 
theory devoted to the analysis of the distribution of control within a network and how the 
behaviour of the system relates to the properties of the components. With few exceptions the 
MCA literature treats networks in steady state. 

MCA has had great success in describing the control and regulation of systems at steady 
state. However, in the analysis of an increasing domain of examples it becomes necessary 
to consider sensitivities along non-steady trajectories. In many systems it is the transient 
or oscillatory behaviour which is of primary interest (e.g. in signal transduction or cell cycle 
regulation). Some extensions of MCA to non-steady behaviour have appeared in the literature 
(|l], [?], O, 12, [ll|]). These papers generalize the steady state results of MCA to special cases 



of dynamic behaviour (e.g. periodic behaviour or trajectories near a stable equilibrium). The 
analysis presented in this paper extends these results by measuring time varying sensitivities 
along arbitrary trajectories. 
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A standard MCA analysis might treat, for example, the effect of increasing an enzyme 
concentration on the steady state value of some related metabolite. Using the analysis described 
below, one can determine the sensitivity to the perturbation throughout the time evolution of 
the system, regardless of the nature of the trajectory. Based on these time- varying sensitivities, 
the Summation and Connectivity Theorems of MCA are extended to give conditions which 
hold for all time. When applied at steady state, these results reduce to the standard analysis 
of MCA. 

The statements in this paper apply to arbitrary system dynamics, including transients, 
convergence to steady state, and oscillations. However, as pointed out in previous papers 
(0, [II]]), sensitivity coefficients derived for oscillating systems must be interpreted with care. 
A discussion on applications to oscillating systems appears below. 

The outline of the paper is as follows. Time-varying sensitivity coefficients are defined and 
derived, and their computation is discussed briefly. The main results of MCA - the Summation 
and Connectivity Theorems - are then extended to the time-varying case. Finally, the sensitivity 
analysis is illustrated by treating some examples, including a model of phototransduction. 

2 Preliminaries 

A network consisting of n chemical species involved in m reactions is modelled. The n- vector 
s is composed of the concentrations of each species. The constant r- vector p is composed of 
the (external) parameters of interest in the model. The m- vector valued function v = v(s, p, t) 
describes the rate of each reaction as a (possibly time- varying) function of species concentrations 
and parameter values. Finally the n by m stoichiometry matrix N describes the network: 
component N, is equal to the net number of individuals of species i produced or consumed in 
reaction j. The network can then be modelled by the ordinary differential equation 

^s(t) = Nv(s(i),p,i) for all t > 0. (1) 

We allow the vector p to contain the initial states of any of the species concentrations as well 
as any other parameters which will directly affect the rates of the reactions (e.g. concentrations 
of enzymes and external effectors). 

We assume that the function v(s, p, t) is continuous and is continuously differentiable in s 
and p for each fixed t. We further assume that for each initial condition s(0) = sq and each 
choice of parameters p the unique solution of ([l]) is defined for all t > 0. We denote this solution 
by s(t, Sq,p), and will use simply s(t, p) or s(t) when no confusion will arise. 

Before embarking on an analysis of system (jl]) it is prudent to first consider any linear de- 
pendencies inherent in the state variables of the system (which will simplify both the analysis 
and the computation). Each conserved moiety in the network corresponds to a linearly depen- 
dent row in the stoichiometry matrix N. We follow the procedure and terminology described 
by Reder [Il4| (see also Q) in the reduction of the system. 

Let no denote the row rank of N. If no = n, then no reduction is necessary. Otherwise, we 
begin by re-ordering the rows of N so that the first uq rows are linearly independent. Let Nr 
be the reduced stoichiometry matrix which results from truncating the last n — uq rows of N. 
Since the truncated rows can be formed by linear combination of the rows of Nr, the matrix 
N can be written as the product 

N = LN R 
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where the n x uq matrix L, called the link matrix, has the form 



A n 

L 



(Here and below the notation I q will be used for the q x q identity matrix). 

The advantage of this decomposition can now be realized. Since each conserved moiety 
allows one species concentration to be determined as a function of the others, we may decompose 
the species vector s into independent and dependent species vectors Sj and respectively. 
Ordering the components of s to match the rows of N, we write s = (sj,Srf) where Sj is an 
no-tuple and is an (n — reo)-tuple. (This involves a minor abuse of notation since s, Sj and 
Sa are all column vectors.) Then equation (Q) can be written as 



d_ 



Si(t) 

s d (t) 



LN R v(s(t),p,t) 



Lo 



Hence 



d 



Sd(t) 



Lo|*(i) 



N R v(s(t),p,t) for all t > 0. 



for all t > 0, 



and so s^(t) — LoSj(t) is an integral of motion of (||), i.e. this difference is constant throughout 
the evolution of the system. We introduce the constant (n — no)-vector T to quantify this 
relationship. Any trajectory of the system satisfies 

s d (t) = L Sj(i) + T for all t > 0, (2) 

where T is defined in terms of the initial conditions as 

T = s d (0) - L Si(0). 

In analysis and computation, attention can be restricted to the independent species Sj, since 
the corresponding results incorporating the dependent vector s d are arrived at immediately 
through the relationship (0). That is, one need only consider the reduced system 



j t Si(t) = N H v(s(t), p, t) = N R v((si(t), L Si(t) + T), p, t) 



for all t > 0. 



(3) 



3 Sensitivity Analysis 

We now present a general sensitivity analysis of the system (|l]). 



3.1 Definitions 



Following [18 1 and O, we make the basic definition. 



Definition 3.1 Given an initial condition s(0) = so and a set of parameter values po, we define 
the time-varying concentration sensitivity coefficients (or concentration response coefficients) as 
the elements of the n x r matrix function R s (-) given by 



R-(t) :=^| P = P0 = Hm 
dp H H Ap^O 



s(t,p + Ap) - s(t,p ) 
Ap 



for all t > 0. 



□ 
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Remark 3.2 These time varying response coefficients can be interpreted in exact analogy to 
the standard (steady state) response coefficients of MCA. The difference is that the response 
to a perturbation along an entire trajectory in now considered, rather than at a particular 
equilibrium state. 

Take for example a system which has an asymptotically stable equilibrium (say s ss ), and 
consider a set of initial conditions (say Sq) which do not match the steady state solution. Letting 
the system evolve from so, we may observe some initial transient, followed by convergence to s ss 
as time tends to infinity. Alternatively, if we make a small perturbation to a system parameter, 
we may observe a different transient, followed by convergence to a different steady state. The 
response coefficient defined above provides a measure of the difference between this "perturbed 
trajectory" and the "nominal" (unperturbed) trajectory at each time t. As time tends to 
infinity, each trajectory will converge to its steady state, and so the response coefficient will 
converge to the steady state response of MCA. 

This time-history analysis is particularly useful when studying systems in which transient 
behaviour plays a key role in the mechanism. For example, in systems performing phototrans- 
duction or action potential transfer, the role of the system is to produce a (transient) spike in 
the concentration of a certain species. The analysis presented here will allow elucidation of the 
effect of parameter perturbations on the behaviour of such a system. 

It should be noted that the response coefficient R s (-) is defined with respect to a particular 
parameter choice po and a particular initial condition so- Should one be interested in comparing 
the system response to perturbations at different times, a separate response coefficient must be 
defined for each such time, since each choice will yield a distinct "initial" state (by re-setting 
time to zero at the point chosen). □ 



Remark 3.3 Since the vector p is allowed to contain the initial states as components, sensi- 
tivity coefficients with respect to initial conditions are included in the above definition. When 
pursuing a steady state analysis, it is often more convenient to consider the total amount of 
each conserved moiety as a parameter rather than the initial conditions of particular species 
(see e.g. ||). However, in this dynamic analysis, perturbation of the initial concentrations of 
different species will have distinct effects on the time-history; it is only at steady state that these 
effects can be described equivalently in terms of perturbations in pools of conserved moieties. 
□ 

In addition to the sensitivities of the species concentrations, it is also of interest to consider 
the sensitivities of the various rates v(s,p, t). In expressing these responses, we conform to 
the MCA usage of the term flux for the rate of a system reaction. (In general, these rates 
do not represent fluxes in the usual sense of the word as a rate of transfer of some quantity; 
the rate of passage of mass through a reaction pathway must also take the stoichiometry into 
consideration.) 

Definition 3.4 Given an initial condition s(0) = so and set of parameter values po, we define 
the time- varying flux sensitivity coefficients (or flux response coefficients) as the elements of the 
m x r matrix function R J (-) given by 

j _ dv(s(t,p),p,t) , _ flv(s,p,t)as(t,p) | dv(s,p,t) 

( ' dp lp=po " ds dp + dp 

= ^ (S ' P 'V (*) + av(S ' P ' t} for all t > 0, 
as op 
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where the derivatives are evaluated at p = po and s = s(f, po). 



□ 



Remark 3.5 The flux response coefficients can be interpreted analogously to the concentration 
response coefficients: R J (i) gives the response in the rates at time t to a perturbation at time 
zero. At steady state, these coefficients reduce to their standard MCA counterparts. □ 



3.2 Computation 



The sensitivity coefficients are defined by a first-order linear ordinary differential equation 
which follows from taking the derivative of (|l|) with respect to p. Dropping some functional 
dependencies for ease of legibility, we have: 



d ds(t) 



N 



fdv(t)ds(t) dv(t)\ 



+ 



dp J 



dp dt \ ds dp 

which becomes, after switching the order of differentiation, 



for all t > 



d ds(t) 
dt dp 



N / av(*)as(t) | dv{t)\ 

\ ds dp dp J 



for all t > 0. 



(4) 



Recall that the sensitivity coefficients are defined with respect to a particular choice of sq = s(0) 
and po- This choice fixes a trajectory s(t) = s(i,so,Po), and hence determines the function 
v(t) = v(s(f),p, t) as well. The initial conditions for (||) (i.e. the components of the matrix 
1^(0)) are immediate from the definition of the sensitivity coefficients: the j, k-th entry of Jjj(0) 
will be zero unless the k-th parameter is the initial condition of the j-th species, in which case 
the initial value of the sensitivity coefficient will be one. Equation (Q) can be solved for the 
concentration sensitivities J-^(-) = R s (')) an d the flux sensitivities can be computed according 
to definition |3.4| . However, we can reduce the burden of computation by taking advantage of 
the dependencies described in equation (|2|). 

Differentiating (g) with respect to p, we find 



d dsj(t) 
dt dp 



N R 
Nr 
Nr 



av(t)^ ( t) + ^^ (LoSj(t)+T) + ^ 

+ dv(t) dT dv(t)~ 



dsi 



dp ds d dp 
( dv(t) + dv(t) L \ ds t (t) 
V dsi ds d °J dp 

dv(t) T dsi(t) dv(t)dT a 

^ — ~ 1 ^ 7^ r 



ds dp dsd dp 
The initial conditions are determined as before. 



ds d dp 

dv(ty 

dp 



dp 

for all t > 0. 



(5) 



Remark 3.6 While there are special cases in which the response coefficients can be derived ex- 
plicitly (as discussed below), in most cases the solution to equation (|j) (or more efficiently (||)) 
must be computed numerically. In performing such a calculation, one must keep in mind that 
the right-hand-side depends on the time- varying value of the species concentration vector s(t). 
A simple computational strategy is to solve (0) and (JsT) as a single set of differential equa- 
tions, determining Sj(-) and simultaneously. This would allow, for example, an adaptive 
integrator to reduce the step-size during computation should either equation demand it. 
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The reader should note that the pair of ODE's are coupled, but that the coupling is in 
cascade - equation (||) is independent of 9s q^ ■ Thus the increase in complexity comes simply 
from having two equations to solve, not from any intertwining of the solutions. 

An algorithm for computation of the time-varying response coefficients has been imple- 
mented in MATL AB and in the biochemical simulator Jarnac . Scripts of the implementa- 
tion are available online ( [http:/ /www.sys-bio.org| ), □ 



3.3 Analytical Solution 

The solution to equation (|5|) can be expressed in terms of the variation of parameters formula 
(derived in standard texts on ODE's, e.g. Q): 

^t) = *(*)*« (0) + Ht) f *-i(r)N R ( + dr for all t > 0, (6) 

dp dp J V ds d dp dp J 

where the matrix- valued function $(•) (called the fundamental matrix for equation (Eh) satisfies 
the homogeneous part of (|5|), i.e. $(■) is defined as the solution of the following initial value 
problem 



N R ^L 

as 



m, m=K- (7) 



Note that @ does not in general provide an explicit formula for the solution to (|B|) (since (|7|) 
does not admit an explicit solution, in general). Nevertheless, knowing the solution takes this 
form can simplify both analysis and computation. In addition, as will be demonstrated in 
Section |4], there are special cases in which (|6|) does indeed provide an explicit solution. 

In the subsequent analysis, it will be convenient to partition the parameter vector p into the 
set of parameters which affect the rates (denoted p„) and the set of initial conditions (denoted 
p s ). From the fact that J^- = J^-(0) = Jj-p = 0, we see that formula @ can be partitioned as 

dsi , , ^ , , ds, , , ^ . , /"* ^ 1 . _ T d\(r) dT , 
— ^ = $(t)_L(0) + *(t / r^rNR-ii-dr 
<9p s <9p s Jo ds d dp s 

t 



^L(t) = / r^rjNR^dr for all t>0. (8) 



To provide a more elegant exposition in what follows, we consider the extension of (|8|) to 
the entire response coefficient matrix R s (-) = fM')- Since s^(t) = LoSj(i) + T for all t > and 
J£ = 0, © gives 

9s, 



-(t) = L $(t) / $ _1 (r)N R ^^(iT for all i > 0. 
« jo dp v 



(9) 



<9p 

Together, (||) and (||) give 

■^-M = L*(*) / rV^Na^dr for all f > 0. (10) 
#Pu Jo <?Pu 

As for the flux-response coefficients, 

£l(t) = ^2L*(t) T ^-WNa^l rfr + 5® for all t > 0. 
<9pu as Jo dp v dp v 
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3.4 Scaled Coefficients 



Computations in MCA are typically done using "scaled" coefficients: derivatives are taken with 
respect to the logarithms of the variables, giving relative (rather than absolute) responses. 
These scaled sensitivities have the advantage that they are dimensionless. Further, they allow 
direct comparison of responses at different states or across different parameters. 

The scaled responses are related to their unsealed counterparts through multiplication by a 
ratio of the values involved. These relationships can be elegantly described at the matrix level. 
Following ||, we define the diagonal matrices 

D s (-) := diag s(-) D v (0 := diag v(s(0, P ,0 D p := diag p 

whose entries are given by the coefficients of the corresponding vectors. The scaled response 
matrices, denoted with a "tilde" (~) are then given by 

dlns(-) 



R s (.) = (D 8 (0) _1 D p R B 



dlnp 



R v (0 = (D v (0)- 1 DPR J (0 = ^^ 



4 MCA Theorems 



We next consider generalizations of the main results of Metabolic Control Analysis to the case 
of time-varying sensitivities. These results (the Summation and Connectivity Theorems, and 
the resulting Control-Matrix Equation) are usually stated in terms of control coefficients and 
elasticities, into which the response coefficient matrices are factored ( fL4| , ||] ) . The time varying 
responses defined above do not, in general, allow such a factorization. However, the responses 
can be interpreted as the result of applying an appropriately defined control operator to an 
elasticity. In this sense, the standard MCA results can be extended to the time-varying case. 
Moreover, we will show that under certain assumptions on the derivatives of the rates, one can 
recover direct generalizations in terms of matrix factorizations (as was done in Rj). 



4.1 Control Operators 

Given a trajectory of system @), we define the substrate- and parameter- elasticities of the 
system as 

These vectors describe component (or "local") sensitivities of the isolated reactions associated 
with the system. 

The steady state Summation and Connectivity Theorems admit elegant mathematical de- 
scriptions due to the fact that the steady state concentration and flux responses can each be 
decomposed into a product of two terms: a matrix of control coefficients and a vector of param- 



eter elasticities. It is clear from equation (10) that the general time- varying response cannot 



be expressed as a product involving the parameter elasticities (since this time-varying term 
appears inside the integral). However, the relation between response and elasticity can be ex- 
pressed as the action of a control operator on the elasticity vector, as follows. (An operator, in 
this case, is an object that maps functions to functions. An introduction to operators appears 
in the appendix.) 
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Definition 4.1 The concentration control operator C s maps m- vector valued functions to n- 
vector valued functions as follows. Given a continuous m-vector valued function y(-) defined 
on the non-negative reals, C s (-)(y(-)) is an n- vector valued function defined by 

C s (t)(y(-)) := L$(i) f $ _1 (r)NR y(r) dr for all t>0. 
Jo 

□ 

In a like manner we define the flux control operator. 

Definition 4.2 The flux control operator C J maps m-vector valued functions to m-vector 
valued functions as follows. Given a continuous m-vector valued function y(-) defined on the 
non- negative reals, C J (-)(y(-)) is an m-vector valued function defined by 

C 3 (t)(y(-)) := e s (t)C s (t)(y(-)) + y(t) for all t > 0. 

□ 

The elasticities and control operators can be scaled in the same manner as the response co- 
efficients. Each of the results described below has an analogous statement in terms of scaled 
quantities, determined simply by including the scaling factors. 

With these definitions in hand, the generalizations of the partitioned response properties 
( |10| , [H ) follow immediately. If the parameter vector is such that the parameters only affect the 
rates (i.e. p = p v ), then 

R s (t) = C s (t)(e p (-)) R J (t) = C J (i)(e p (-)) for all t > 0. 

We next indicate how the Summation and Connectivity Theorems can be interpreted in the 
light of these definitions. We will also show how the action of these operators reduces to matrix 
multiplication when the system is at steady state. 

4.2 The Summation Theorem 

The Summation Theorem can be stated in terms of the control operators as follows. 

Theorem 1 (Summation Theorem) If the m-vector k lies in the nullspace o/N (i.e. Nk = 
0, and so NR,k = as well), then (using the simplified notation described in the appendix) 

C s (i)(k) = and C J (t)(k)=k for all t > 0. 
Proof. The result follows from the definitions of the control operators. If NrIc = 0, then 

C s (t)(k) = L$(i) f *- 1 (r)N R kdr = for all t > 0, 
J o 

and 

C J (t)(k) = e s C s (t)(k) + k = k for all t > 0. 

I 

To state the result in a more standard form, we allow the control operators to act on matrices 
as follows. If ri(-), . . . , r n (-) are m-vector valued functions, then we interpret 

C s (t)([ ri (.) . . . r n (-)D := [C s (t)( ri (.)) . . . C%t)(r n (.))] . 

With this notation, a direct corollary of the Summation Theorem is the following. 
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Corollary 4.3 // the matrix K has columns which form a basis for the nullspace o/N, then 
C s (t)(K)=0 and C J (t)(K) = K for all t > 0. 



The Summation Theorem is normally stated as a property of the control coefficients of a 
system at steady state. However, the steady state version can be given an equivalent formulation 
which does not make reference to control coefficients. 

Proposition 4.4 (Summation Theorem — steady state) Suppose the system is at steady 
state, so we may drop time dependencies on all terms. If the parameter vector is such that the 
parameters only affect the rates (i.e. p = p„) and each column of the matrix e p = |^ lies in 
the nullspace of N, then 

R s = and R J = e p . 

■ 

Stated in this way, the Theorem can be extended directly to time-varying responses. From 
Theorem [|, we have the following. 

Corollary 4.5 (Summation Theorem — alternative statement) If the parameter vector 
is such that the parameters only affect the rates (i.e. p = p v ) and each column of the matrix 
£ p (t) = jjjj(£) lies in the nullspace o/Nr for each t > 0, then 

K s (t) = and R J (t) = e p (t) for all t > 0. 



4.2.1 Special Case: e p constant 

In special cases, Theorem |] reduces to more familiar variants. A simplification which can be 
made to the general case described above is to assume that some derivatives of the rate law 
v(-, •, •) are constant throughout the evolution of the system. This was the case considered by 
Heinrich and Reder in 0, where it was argued that for systems near a stable equilibrium, it 
may be reasonable to approximate the derivatives of v by their values at the equilibrium. 

If one restricts to the special case where the parameter elasticity e p is constant throughout 
the evolution of the system, the control operators need not be defined to act on arbitrary 
functions, but only on constant m-vectors. The action of the operators can in this case be 
described simply as matrix multiplication, e.g. for any constant m-vector y, 



C s (t)(y) 



L$(t) f $ -1 (r)N R dr 
Jo 



for all t > 0. 



In this case, through a slight abuse of notation, time varying control coefficients can be 
defined by 



C s (t) := L$(t) /V^NRdT 
Jo 

C 3 (t) := e 8 (i)L*(i) f $ _1 (r)N R dr + I r 

Jo 



9 



Here, we can state a more standard Summation Theorem: if the matrix K has columns 
which form a basis for the null space of N (equivalently of Nr.), then for all t > 0, 

C s (t)K = 
C J (i)K = K, 

where the left-hand-side is interpreted as a matrix product. 
4.3 Connectivity Theorem 

The Connectivity Theorem for control operators is as follows. 

Theorem 2 (Connectivity Theorem) Applying the control operators to the function £ s (-)Jj = 

dv i 
ds [ 



)L yields, for all t > 0, 



C s (t)(e s (.)L) = L($(t)-I no ) 
C 3 (t)(e s (-)L) = E 8 (t)L*(i). 

While the right-hand-sides of the equations above may not look familiar, the reader should note 
that if <£(t) = then these reduce to more standard connectivity relations. In what follows we 
will consider a case where <J>(i) approaches zero as steady state is reached, so that the classical 
MCA result is recovered. 

Proof. The definition of the concentration control operator gives 



C s (t)(^(-)L) = L$(i) jf r 1 ^ ^^ LdT for a11 * ^ °- 

From equation (|7|), we have that 

" d 



dt 



*(t) 



as 



for all t > 0. 



(11) 



(12) 



Moreover, from the fact that 
d, d 



= H l '" = J^WW 



* _i (t)+#(t) 



we see that 



AcD-i(t) = _*-i( t ) 



d ^ , « 
d^ } 



for all i > 0. 



Thus, from (|T|) and 

C s (t)(^(-)L) = L*(t) j\-i( T )N n ?^p-LdT 

= L$(i) /% _1 (t) F-^-$(t)1 * _1 (T)dT 
Jo L« r 

= -L$(i) / — $ _1 (r)dr for all t > 0. 
Jo « T 
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Application of the Fundamental Theorem of Calculus gives 

C s (t)(^(-)L) = — L$(t)[$ -1 (t) — <& -1 (0)] 

= L($(i) -I no ) for all t > 0, 

as <]? -1 (0) = I no = $(0). With the definition of the flux control operator, we conclude 

C J (t)(e s (-)L) = e 8 (*)L(*(t)-I no ) + E,(t)L 
= e s (t)L$(t) for all t > 0. 

I 

The Connectivity Theorem is typically stated as a property of the steady state control 
coefficients. As is the case for the Summation Theorem, an equivalent formulation can be given 
which does not make reference to control coefficients, as follows. 

Proposition 4.6 (Connectivity Theorem — steady state) Suppose the system is at steady 
state, so we may drop time dependencies on all terms. If the parameter vector is such that the 
parameters only affect the rates (i.e. p = p„) and the elasticities satisfy e p = e s L, then 

R s = -L and R J = 0. 

■ 

There is a direct extension to the general case. From Theorem ||, we have the following. 

Corollary 4.7 (Connectivity Theorem — alternative statement) If the parameter vector 
is such that the parameters only affect the rates (i.e. p = p v ) and the elasticities satisfy £ p (t) = 
e s (i)L for each t > 0, then 

R s (t) = L($(t) -I no ) and R J (t) = e s (t)L*(t) for all t > 0. 



4.3.1 Special Case: e s and e p constant 



Again, in the special case where the derivatives of v(-, •, ■) are constant throughout the trajec- 
tory, the analysis simplifies. 

If e s = 4j£ is constant throughout the evolution of the system, the fundamental matrix 
defined by equation (^) can be described explicitly as the matrix exponential 



*(*) 



3 N R£s Lt 



for all t > 0. 



In this case (|6|) provides an explicit formula for the sensitivity coefficients. 

Heinrich and Reder considered the case in which both e p and e s are constant along the 
trajectory of the system in (tJ. In that case, the control operators again need only be defined 
on constant vectors, their action reducing to matrix multiplication. That is, given any constant 
m-vector y, 



C s (t)(y) 



Le NR£sL * £ e- NR£sLT N R air 
. Jo 



for all t > 0. 



(13) 
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The matrix N R e s L is the Jacobian of equation (0). In the case where it is nonsingular, ( |l3| ) 
further reduces to 

C s (t)(y) = [L( e N ^ L '-I no )(N R e s L)- 1 N R ]y for all t > 0. 

This assumption of nonsingularity is rather standard (see, e.g. |jT3f ); in particular, it holds when 
the analysis is being carried out near an asymptotically stable equilibrium point. 

Under these assumptions the control coefficients take the form of the time- varying matrices 

C s (t) := L(e N ^ Lt -I no )(N R£s L)- 1 N R 

C 3 (t) := es L(e NR£sL '-I„ )(N R£s L)- 1 N R + I m , 

and the Connectivity Theorem reduces to the form presented in j?]]: for all t > 0, 

C s (t)e s L = L(e N ^ L '-I no ) 
C 3 (t)e s L = e s Le NR£sL *. 

If such a a trajectory further satisfies the condition that it is approaching an asymptotically 
stable steady state, then it will be the case that 

lim e NR£sL * = 0. 

Then, taking limits, we find that the control coefficients reduce to their standard forms as 
the steady state is approached, 

limC s (t) = -L(N R e s L)- 1 N R 

t— »oo 

limC J (t) = - es L(N R£s L)- 1 N R + I m . 

The time-varying Summation and Connectivity Theorems described above reduce to the stan- 
dard statements in this case, as described in Jq|. 



4.4 Control Matrix Equation 

Together, the steady state Summation and Connectivity Theorems provide a relationship be- 
tween the elasticities and control coefficients. An elegant description of this relation is the 
Control Matrix Equation || . This equation has a direct generalization in terms of the control 
operators described above. 



Corollary 4.8 If K is a matrix whose columns form a basis for the nullspace of N, then, 
making use of the notation for vectors of operators introduced in the appendix, 



c J (t) 

C s (t) 



( K -e s L ) 



K 





-e.(t)L*(t) 

L(In " *(*)) 



for all t > 0. 



□ 
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Figure 1: Pathway 
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Figure 2: Perturbation in k\ 

5 Applications and a Caveat for Oscillating Systems 

We present applications of the above analysis to three models. The first is a "toy" model 
- a stable two-species pathway. This analysis provides a straightforward illustration of the 
interpretation of time-varying response coefficients. The second example is a simple oscillating 
system. In light of this example we provide a brief discussion of the interpretation of sensitivity 
coefficients of oscillatory systems. Finally, we provide a "realistic" example - an analysis of a 
phototransduction model based on work of Reike and Baylor Q. Here the power of a time- 
varying sensitivity analysis is demonstrated, since we are able to make statements about the 
effect of parameter variations on the critical transient behaviour displayed by the model. 

5.1 Basic Pathway 

Consider the simple pathway shown in Figure |l]. Assume mass-action kinetics for each of the 
reactions. We will model the system with the following parameter values: vq = 4, k\ = 3, 
= 0.5, and k<i = 2. Since the flow into Si is fixed, the system is asymptotically stable to 
the equilibrium (Si, S2) = (| , 2). We will consider the unsealed response coefficients defined for 
initial conditions at the steady state (recall R s (-) depends on the nominal parameters values and 
the initial state). Thus in this case the nominal trajectory will be the steady state trajectory 
(Si(t),S 2 (t)) = (§,2) for all t > 0. 

First consider the sensitivity to a perturbation in k\ (shown in Figure §). The steady state 
response of S2 to such a perturbation is nil (since at steady state, S2 = j^). However, the 
analysis shows (and intuition concurs), that an increase in k\ will produce a transient increase 



13 



1 



I — Sensilivily - S1 I 
| Sensitivity - S2 | 



0.5 1 1.5 2 2.5 3 3.5 4 

Time 



Figure 3: Perturbation in k-i 
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Figure 4: Perturbation in Si(0) 



in 52, which fades gradually as Si decreases to accommodates the perturbed parameter. The 
response in S± grows to reach its steady state value - a small overshoot can be seen at about 
time t = 1. Not surprisingly, the situation for perturbations in is similar (but reversed), 
with a transient decrease in S2 and a positive steady state response in S\ (Figure [|). One can 
read from the graphs how the effect on each concentration changes over time. For instance, the 
effect of these changes on S2 is felt most strongly at about t = 0.4. 

Finally, we consider the effect of a perturbation in the initial value of S\ (Figure Again, 
the results are as expected, with transient effects of the perturbations "washing through" the 
pathway - the result being no effect on the steady state solution. 

Having considered the effect of perturbations around a steady state, we now turn our at- 
tention to a system whose time-behaviour is more complex. 

5.2 Oscillatory Systems 

We consider another simple two species pathway, based on a model which has been employed 
in investigations of glycolytic oscillations || . The network is shown in Figure ||. The rates are 
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Figure 5: Oscillatory System 




Figure 6: Perturbation in S'i(O) 



calculated by mass-action, with the activation by S2 appearing as a multiplicative factor (i.e. 
the rate of production of S2 is given by fei5i(l + Sf))- For parameters in a certain range, the 
concentrations of Si and S2 follow oscillatory trajectories; the positive feedback from S2 drives 
the oscillation. 

We choose nominal parameter values of vq = 8, fei = 1, &2 = 5, k$ = 1, and q = 3, 
and nominal initial values of (S'i(O), ^(O)) = (1,3). The oscillatory behaviour is a limit cycle 
(cf. H), and so any initial conditions in a neighbourhood of these points will yield similar 
trajectories. 

We first consider perturbations in initial values. The limit cycle behaviour is stable; after 
a perturbation of the initial conditions, the system will tend to the limit cycle as time tends 
to infinity. However, unlike the previous case where convergence to the same behaviour meant 
that the response shrank to zero as time moved on, for this model (and for systems with 
limit cycles in general), the phase of the oscillations depends on the initial values. Thus the 
perturbed trajectory is out of phase with the nominal trajectory. The result is a periodic 
response - oscillating with the same period as the limit cycle and indicating the difference 
in the trajectories at each point in time. The response for perturbations in S'i(O) is shown 
in Figure H along with the nominal trajectory. The response for perturbations in 62(0) (not 
shown) is similar. 

Care must be taken in extending this set of ideas to perturbations of kinetic parameters of 
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an oscillatory system. In general, the results of such analysis may not be useful. The response 
coefficients for perturbations in vq are shown (again, along with the nominal trajectories) in 
Figure [7[ While the sensitivities match the period of the system in a quasi-periodic behaviour, 
they oscillate more and more widely, and continue to do so as time tends to infinity. As a result, 
these responses are not especially useful when trying to predict the effect of finite perturbations 
on a real system. The trouble, as pointed out in |3], pT| , is that after a perturbation in a kinetic 
parameter, the system again exhibits a limit cycle behaviour, but this time with a different 
period. When two trajectories follow a similar cycle with different periods, they will always 
reach a point when they are completely out of phase, no matter how small the difference in 
their periods. As the perturbations get smaller, the difference in period shrinks as well, so that 
it takes longer and longer for this "point of maximum difference" to be reached. However, as 
long as the perturbation is nonzero, that point will always be achieved. Since the response 
coefficient is defined as the limit of the ratio of the difference in trajectories to the size of 
the perturbation as the perturbation goes to zero, we find that as time tends to infinity, the 
response approaches this "maximum difference" divided by zero, that is it diverges. 

In UH and p| , some clever strategies have been devised for interpreting the sensitivities of 
oscillatory systems. Kholodenko et al. treat the case of asymptotically stable systems under 



the influence of periodic external forces in [11]. Such systems exhibit limit cycles whose period 
is set by the external force. In this case, perturbations in kinetic parameters do not result in 
changes in the period, and so the response coefficients are more meaningful. In Q, Demin et 
al. treat perturbations of autonomously oscillating systems by measuring the response in the 
Fourier coefficients of the trajectory. Fourier control coefficients are introduced, which give a 
useful interpretation of sensitivities for any periodic systems. 

Having considered the interpretation of response coefficients for steady state and oscillating 
systems, we now treat the most interesting case - a system in which the behaviour of interest 
occurs during transients. 

5.3 Phototransduction 

Phototransduction is the process by which organisms convert light signals into nerve signals. As 
in all signal transduction pathways, the steady state behaviour of the system is less interesting 
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Figure 8: Phototransduction Pathway 



than its transient response. We will illustrate a sensitivity analysis of the single photon response 



in a simple model of phototransduction loosely based on previous modeling work in ||] and [13]. 

The response mechanism is modelled as follows. The system is at steady state in the 
absence of light input. A photon of light activates a rhodopsin molecule to start the system 
response. Activated rhodopsin (Ra) in turn activates a heterotrimeric G protein (G). The 
activated G protein a-subunit (Ga) dissociates from its beta-gamma subunits (P"f) and then 
binds with a cGMP phosphodiesterase (PDE) to form the active PDEa-Ga complex. The PDE 
enzyme catalyzes the degradation of cGMP in the cell. Activation by Ga greatly enhances the 
catalytic activity of PDE, resulting in a reduction in cellular levels of cGMP. The cGMP-gated 
ion channels in the cell membrane then close, causing a net efflux of calcium ion (Ca), a net 
outward current (J), and a corresponding hyperpolarization of the cell's membrane. 

The deactivation phase of the response is modelled as follows. Activated G protein a 
subunit, either free (Ga) or bound (PDEa-Ga), is converted to the deactivated form (Gd). 
Finally, the deactivated G protein (Gd) recombines with the beta-gamma subunit (/3j) to form 
the heterotrimeric G protein (G). 

The model incorporates four independent variables (G, Ga, PDEa-Ga, and cGMP) along 
with five dependent variables (Gd, PDE, /?7, Ca, and J) and several parameters as indicated 
below. Note, J denotes not the actual current but rather the magnitude of the deviation from 
the resting current (DarkJ). The input to the system is the time- varying level of Ra. The 
network is shown in Figure ||. The dynamics are given by 

vi = k Ga Ra 
V2 = k Gd Ga 
^3 = kai /?7 Gd 

«4 = kpoEi Ga PDE — kpDEim PDEa-Ga 
v§ = kpDE2 PDEa-Ga 
k c G 



1+ i Ca 



V7 = k cGd0 cGMP + k cGd cGMP PDEa-Ga 
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For simulations, we choose nominal parameter values and initial conditions as follows 



1 2 

/im 



kca = 1-0 ms 1 
feci = 0.002 ms" 1 molecules 
kpDEim = 0.000001 ms" 1 
k cG = 0.058 mM ms _1 
q = 2.0 

k c Gd = 0.000004 ms _1 molecules 
k Cl = 0.72 

Gtot = 3000 molecules /im -2 
P^tot = 3000 molecules fimT 2 

G (0) = 3000 molecules jum" 2 
PDEa-Ga (0) = molecules jimT 2 



k Gd = 0.005 ms" 1 
kpDEi = 0.0003 ms" 
kpDE2 = 0.005 ms~ 
K g = 0.2 mM 
kcGdo = 0.002 ms" 1 
//m 2 kc a = 1.44 

#0.5 = 4 mM 



1 molecules 1 /jm 2 



PDEtot = 500 molecules fim 2 
DarkJ = 12.0 pA 

Ga (0) = molecules //m -2 
cGMP (0) = 4.0 mM 



The single activated rhodopsin molecule excited by a single photon is represented by a 
square 100ms pulse, 

( < t < 100 
Ra (t) = < 1 100 < t < 200 , 
( 200 < t 

measured in molecules /im" 2 . 

We first consider the effect of a perturbation in the parameter k Ga , which describes the 
effect of activated rhodopsin on the activation of G. The nominal trajectories for G, Ga and 
PDEa-Ga are shown in Figure |9] and Figure 10. (There is no activity until time t = 100 since 



the system is at steady state while the input is zero). Figure 11 shows the absolute (unsealed) 
sensitivity of G, Ga and PDEa-Ga to perturbations in k Ga - As expected, an increase in kc a 
causes a decrease in the level of G and an increase in the levels of Ga and PDEa-Ga. These 
responses track the trajectories themselves, with the largest response being observed when each 
species is farthest from its "resting" (steady) state. 

The primary signal of interest in this system is the level of current (J) produced by the 
input. The time history of the current produced by the nominal input and parameter values 
is shown in Figure 12. The effect of perturbations in three different parameters on J are 
shown in Figure [l3|. Since this analysis is meant to compare the effects of the changes, the 
relative (i.e. scaled) responses are shown. The relative strengths of perturbations in the three 
parameters k Ga , k G d and kpDEi are immediate. For instance, to increase the response in J, it 
is clear that an increase in k Ga will have far more impact than an equal (relative) increase in 
kpDEi or decrease in k G d- 
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Figure 9: Nominal trajectory for G 




Figure 10: Nominal trajectory for Ga, PDEa-Ga 
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Figure 13: Relative sensitivities of J 
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6 Conclusion 



Since its introduction, Metabolic Control Analysis has proven to be a useful tool in discovering 
the distribution of control in biochemical systems at steady state. In this paper, the compu- 
tation and interpretation of system sensitivities has been considered for arbitrary trajectories. 
The main results of MCA (Summation and Connectivity Theorems) have been shown to have 
a valid interpretation not just at steady state, but throughout the system dynamics. In analyz- 
ing the distribution of control for systems in which transients are of interest, the time varying 
sensitivity may prove to be a valuable tool in determining system behaviour. 
Acknowledgment The authors would like to thank Tau-Mu Yi for providing the phototrans- 
duction model and for many useful discussions. 



A Appendix: A Primer on Operators 

A function from M into M is defined as a rule which assigns to each real number^ a unique real 
number in its range, for example, the function h(-) defined by h(t) = t 2 for each t € M. This 
can be generalized by allowing functions to act on pairs of numbers, or n-tuples in general (e.g. 
g(x,y,z) = x 2 - yz). 

The concept of a function can be extended to yield an operator, which takes a function as 
its argument. For example, we could consider the operator T(-) defined by integration: 

n/(-))= [ 1 f(s)ds 

Jo 

for each function /(•) defined on M (provided the integral of / exists). Define the functions h 
and r by h(t) = t 2 and r(t) = cos(t). Then 



T(h(-)) 



[ s 2 ds = -, and T(r(-)) = f cos(s) ds = sin(l). 
Jo 3 J 



Another example of an operator is function evaluation. Define G by G{f{-)) = /(5) for each 
function /(•). Then G(h{-)) = 25 and G(r(-)) = cos(5). 

The operators T and G defined above are commonly referred to as functionals, since they 
map functions to numbers. More generally, by allowing an operator to have a second, scalar ar- 
gument, we can construct operators which map functions to functions. For example, extending 
T to depend on a scalar argument t, we could define T by 

%/(•))= f f{s)ds. 
Jo 

An alternative notation is T(t)(/(-)) = T(i, /(•)). Thus, given a real-valued function, T returns 
a function defined for t > 0: 

rt ^3 rt 

T(t)(h(-)) = / s 2 ds = — and T(t)(r(-)) = / cos(s) ds = sin(i). 
Jo 3 J 

As another example, we could modify G by defining G{t)(f(-)) = f(2t). Then G(t)(h(-)) = (2t) 2 
and G(t)(r(-)) = cos(2t), for all t G M. 

"The symbol R denotes the set of real numbers. 
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As described in the main text, we can extend the action of operators to row vectors in a 
component- wise fashion, for example 

T(t)([h(.) r(-)]) = [$ sin(i)]. 

Moreover, we can combine operators column-wise, yielding the notation 



T(t) 
G(t) 



([ MO r(-) ]) 



V sin(t) 
(2t) 2 cos(2i) 



Particular cases of operator actions which appear in the main text are: 

• constant functions as arguments to operators. For example, defining the function k(-) by 
k(t) = 3 for all t € R, we have 

T(*(0) = 3, f (*)(fc(0) = 3t. 

Using a simplified notation, we can write 

T(3) = 3, T(*)(3) = 3t; 

• operators which act by multiplication. For example, define the operator H by H(t)(f(-)) = 
3tf(t). Then 

H(t)(h(-)) = 3t 3 H(t)(r(-)) = 3tcos(t) H(t)(3) = 9t. 
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